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We investigate the evolution of entanglement in the Fenna-Matthew-Olson (FMO) 
complex based on simulations using the scaled hierarchy equation of motion (HEOM) 
approach. We examine the role of multipartite entanglement in the FMO complex 
by direct computation of the convex roof optimization for a number of measures, 
including some that have not been previously evaluated. We also consider the role 
of monogamy of entanglement in these simulations. We utilize the fact that the 
monogamy bounds are saturated in the single exciton subspace. This enables us to 
compute more measures of entanglement exactly and also to validate the evaluation 
of the convex roof. We then use direct computation of the convex roof to evaluate 
measures that are not determined by monogamy. This approach provides a more 
complete account of the entanglement in these systems than has been available to 
date. Our results support the hypothesis that multipartite entanglement is maximum 
primary along the two distinct electronic energy transfer pathways. 



I. INTRODUCTION 

Photosynthesis is one of the most common phenomena in nature. However, the details 
of photosynthetic processes are still under investigation. Recent experimental results show 
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that long lived quantum coherences are present in various photosynthetic complexes [TH3]. 
One such protein complex, the Fenna-Matthews-Olson (FMO) complex from green sulphur 
bacteria [4], has attracted a great deal of experimental and theoretical attention due to 
its intermediate role in energy transport. The FMO complex acts as a molecular wire, 
transferring the excitation energy from the light-harvesting complex (LHC) to the reaction 
center (RC) |1]-|7]. In 2007, Engel et al observed long-lasting quantum beating over a 
time scale of hundreds of femtoseconds by two-dimensional nonlinear spectroscopy. Evidence 
for quantum beating, and therefore long lived quantum coherence, was also found at room 
temperature [9j. 

The transport of electronic excitations through the protein complex of FMO is an exam- 
ple of energy transport in an open quantum system. The oscillations of the nuclear positions 
provide a bath or an environment for the electronic excitations. Since 2007, several theoreti- 
cal frameworks have been developed to model this phenomenon. For example, Aspuru-Guzik 
et al [T014T2] introduced a non-Markov approximation based on the Lindblad formalism to 
investigate the effects on the efficiency of photosynthesis of the combination of quantum co- 
herence and environmental interaction. Meanwhile, Ishizaki and coworkers |T3| [H] utilized 
the hierarchichal equation of motion (HEOM) approach to reproduce the population beating 
in the FMO complex successfully at both cryogenic and physiological temperature. More 
recently, Zhu and coworkers |15] introduced the scaled HEOM approach [16] for studying 
the robustness and quantum coherence in FMO complex. The scaled HEOM approach has 
been proved to provide a reliable simulation result with considerable reduction in compu- 
tational requirements. Using the HEOM equations, Rebentrost and Aspuru-Guzik showed 
that the non-Markovianity of the system is near-maximal at physiological conditions |17] . 
Recently, many other approaches for the numerical computation of the time evolution and 
quantum features of this system have made FMO a target for matchmarking of methods for 
simulating open quantum systems [T8H34] . 

Besides modeling of population and coherence observed in experiment , these models also 
enable computation of the entanglement evolution [3SJ. The first study of entanglement in 
biological excitons was [36], which studied the dynamics of the negativity [371 138] for a pair 
of chromophores coupled to a non-Markovian environment. Subsequent studies considered 
more chromophores, different excitation mechanisms and different entanglement measures. 
We briefly review this work here, for a more complete overview we refer the reader to a 
recent review [39] . In a recent study, Mukamel made a distinction between some apparent 
entanglement effects associated with the linear response which can be eliminated by coor- 



3 



dinates transformation and genuine entanglement that is fundamentally quantum in nature 
|4U| . Recently, Engel et al found a direct evidence of quantum transport in FMO complex 

m- 

In |42] two measures of entanglement relevant to FMO are defined. The first measure 
is the concurrence between chromophore i and chromophore j. The concurrence is a well- 
known measure of entanglement between two two-level systems, and can be computed in 
closed form even for mixed states, and in the case of a density matrix restricted to the single 
exciton subspace takes the simple form = 2\pij\ |^ 115] . The second measure defined 
was a global measure related to the relative entropy of entanglement, defined by; 

N 

E[p] = -J2lnpu-Sip) (1) 

where S{p) = — Trplnp is the von Neumann entropy of the state p. This measure is the 
relative entropy of entanglement specialized to the case where states only have support in 
the zero and one exciton subspace. The definition of the relative entropy of entanglement is 

= min Tr(plnp — p In cr) (2) 

a 

where the minimization is taken over all separable states a. In the case of states restricted 
to zero or one excitons, the set of separable states becomes simply the set of diagonal density 
matrices, and so this minimization can be performed exactly, yielding the expression ([T]). We 
refer the reader to the supplementary materials of [l^ for more details. Both of the measures 
computed in |42| rely on the fact that, in the single exciton subspace, coherence (meaning 
nonzero off diagonal elements of the density matrix in the standard basis) is necessary and 
sufficient for entanglement. Both concurrence, the relative entropy of entanglement and an 
entanglement witness introduced in |42j show this clearly. 

We introduce the notation that the bipartite entanglement between subsystems A and 
B is denoted A\B, and when a subsystem consists of a set of chromophores we indicate 
it by a string of labels, so 12|367 is the entanglement between the subsystem composed of 
chromophores one and two (12), and the subsystem composed of chromophores three, six 
and seven (367). 

The two measures considered in |l2] are computed for an initial excitation at site one 
or six, at both 77K and 300i^, to probe both physiological conditions and the conditions 
of ultrafast spectroscopy experiments. For the system initialized with an exciton at site 1, 
they show the pairwise entanglement 1|2, 1|3, 1|5 and also the pairwise entanglement 3|4. 
Finite entanglement was found between all pairs of chromophores in |42] - over distances 
comparable to the size of the FMO complex - < 30A. 
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The logarithmic negativity is the only measure that is readily computable for all states, 
and in the case of states restricted to the single exciton subspace it may be computed across 
any cut of the set of seven chromophores into two subsets jUHlS]. Caruso et al. computed 
the logarithmic negativity across six cuts 1|234567, 12|34567, 123|4567, 1234|567, 12345|67 
and 123456 1 7 in a simulation in which a single excitation was injected into site one |16]. The 
entanglement of site one with the rest 1| 234567 exhibited the largest peak value, with large 
oscillations taking it below the entanglements across the other cuts. This may be understood 
as the generation of entanglement as the delocalization of the injected exciton across the 
complex. In subsequent work, the logarithmic negativity was also computed (across the same 
cuts) for simulations in which direct injection of a single exciton is replaced by simulation of 
thermal injection and laser excitation. In the case of thermal injection the entanglement is 
reduced by a factor of roughly 50, concomitant with a suppression of coherent oscillations. 
In the case of simulated laser excitation a large pulse of entanglement is observed, lasting 
about 0.15 ps. 

In |47| Fassioli et al move from consideration of the presence of entanglement in models 
of FMO to characterization of its functional role in transport. It is in this context that 
the variety of entanglement studies carried out could connect with functionality and delo- 
calization ideas from physical chemistry. Those authors introduce an entanglement yield, 
based on an entanglement measure which is a sum of the squared concurrences or "tangles" 
(defined below) over all pairs of chromophores. 



Because of monogamy of entanglement their measure is bounded above by a sum of the 
tangles of each chromophore with the rest. 



This upper bound is equal to 2/7 times the Meyer- Wallach measure for the seven chro- 
mophore system j48j . In fact it is known that monogamy bounds are saturated in the single 
exciton manifold [l9] and so the measure Et of Fassioli is in fact exactly equal to 2/7 times 
the Meyer Wallach Measure. Interestingly, those authors point out a connection of this 
measure, and hence of the Meyer- Wallach measure, to a measure commonly used by the 
physical chemistry community of exciton delocalization: the inverse participation ratio |50] . 

To make a connection between entanglement and transport Fassioli et al. |47) define an 
entanglement yield - the integral of the entanglement (as given by a sum of pairwise tangles) 




(3) 




(4) 



n 
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weighted by the probabihty density for exciton absorption by the reaction center. This quan- 
tity is normahzed by the quantum yield: the total probability that the exciton is trapped 
by the reaction center. The contributions to this quantity were divided into donor-donor, 
donor-acceptor contributions, where chromophores 1,2 and 5,6 are designated donors and 
chromophores 3 and 4 are acceptors. This study showed that entanglement peaks on a 
timescale relevant for transport, for simulations in which the initial exciton is localized on 
site one or site six. In particular those authors observe an inverse relationship between en- 
tanglement among donor sites and quantum efficiency, suggesting that entanglement among 
the donor chromophores (1,2 and 5,6) may be tuned to achieve the desired quantum ef- 
ficiency. The authors of [47J also introduce the idea of direct and indirect pathways - an 
indirect pathway involving transfer through chromophore seven. 

In ini] a number of distinct measures of quantum correlation were computed. The quan- 
tum mutual information, quantum discord and single-excitation relative entropy of entangle- 
ment with respect to bipartite cuts 3|16, 12|3 and 3| 124567 were computed. These authors 
extended the work of 0^ by proving a simple formula for the relative entropy of entangle- 
ment across any bipartite cut for states restricted to the single exciton subspace. 

It is the goal of the present work to further investigate the relationship of multipartite 
entanglement to the different transport pathways in the context of the HEOM model pre- 
sented below. The paper is organized as the follows. In Section [IT] the detailed theoretical 
framework of the scaled HEOM approach is introduced. In Section |III| the method used to 



compute the convex roof and hence obtain the entanglement is given. Section IV contains 
simulation results of multi-partite entanglement evolution. We compute cases where the 
entanglement can be determined exactly through the monogamy bound in order to validate 
our convex roof method, and the use the convex-roof optimization to obtain measures that 
are not determined by the convex roof. We close the paper with some conclusions and 
directions for future work. 



II. METHOD: SCALED HIERARCHICAL EQUATIONS OF MOTION (HEOM) 

The structure of the FMO complex was originally analyzed by Fenna and Matthews 
[1]. The FMO complex consists of three identical monomers arranged in a C3 symmetric 
structure. Each monomer is formed from seven bacteriochlorophylla (BChla) molecules. 
These molecules are the "sites" referred to in the rest of the paper. Each monomer works 
independently in the FMO complex. Experimental results show that site 1 and 6 are close 
to the light Harvesting complex (LHC) and site 3 and 4 are next to the reaction center 
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(RC) I1H7]. 

For all models used in the present paper, the Hamiltonian of the FMO complex and its 
interaction with the environment is taken to be: 

"H = Hs + T^B + 'HsB (5) 

N 

N N NjB p2 , 

j=l j=l 5=1 

N N N 

with Vj = and Fj = ^ " 

The terms 7/5, "H^ and "H^b describe the Hamiltonian of the system, the bath, and 
the system-bath coupling respectively. The Hamiltonian is written in the single excitation 
subspace, so that the basis states \j) in Eq. [6] denotes that the j-th site is in its excited 
state and all other sites are in their ground states. The energy of site j is denoted by ej 
and Jjk is the electronic coupling between site j and k. N is the number of sites, so that 

= 7 for the FMO complex. For the thermal bath Tin, the harmonic oscillator model is 
applied. We assume that each site is coupled to the bath independently. The parameters 
rrij^, Uj^, Pj^ and Xj^ are mass, frequency, momentum and position operator of the harmonic 
bath associated with the j-th site respectively. The parameter Cj^ in Eq. |8] represents the 
system-bath coupling constant between the j-th site and ^-th phonon mode. The system 
and bath are assumed to be decoupled at t = 0. 

We can obtain the time evolution of the system density matrix p (t) by tracing out the 
bath degrees of freedom p (t) = TiB[ptotit)] = Tr^ [e'^^*^" Ptot (0) e'^*/'*] . The correlation 
function for a phonon bath can be written as 

00 

1 r -iu)t 
CAt) = - J du.J.iu).^-^ (9) 

—00 

J,H = E^;|^5(^-^.c) (10) 

with /3 = y^sT . We assume that Jj [u) is the same all sites, Jj (u) = J (u) V js. We con- 
sider the time evolution of the system density matrix both with and without environmental 
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interaction. For the isolated system, we set J (w) = and the time evolution of the density 
matrix for the system is given by: 

J^pit) = -{[^s. P{t)] (11) 

One approach to the computation of the time evolution of the system density matrix is the 
hierarchichal equation of motion (HEOM) approach, originally developed by Ishizaki and 
Fleming . We use the scaled HEOM approach for reasons of computational efficiency [T51 



In the scaled HEOM approach, the original spectral density function J (w) (Eq. 10) is 
replaced by a Drude spectral density function J {u) = ^ zj^rp^ where A is the reorganization 
energy and 7 is the Drude decay constant. Then the correlation function in Eq. |9] can be 
expanded as 



C, (t > 0) = 5^ Cfc ■ e- 

*;=0 



with V, 



— 7, which is the Drude decay constant, Vk = ^ when k ^ 1 and Vk is known as 



the Matsuraba frequency. The constants Ck are given by 



Co 



Ck 



V7 
2 



cot 



^2 _y 



for k ^ 1 



Using the scaled approach developed by Shi and coworkers [16j and applying the Ishizaki- 
Tanimura truncating scheme [521 [53] to the density matrix, the scaled density operator 
becomes: 



d 



N K N 

-- [Hs, Pn]-^^ njkVk ■ Pn-i^ 

j=l k=0 j=l 
TV 00 N K 



n 



jk 



+ l)\Ck\ 



'jk 



j=l m=K+l ''■^l 1— n 



j=l k=0 



jk 



jk 



VA (12) 



where the global index n denotes a set of nonnegative integers n = {ni,n2 yn^} = 

{{niQ, rill ■ ■ ■ 5 ^ik} ■ ■ ■ {nm,fT'Ni ■ ■ ■ jIT-nk}}- The symbol refers to a set in which the 
number n^fe is modified to nj^ ± 1 in the global index n. The sum of is called the 
tier (A/"), A/" = ^jj^njk- The global index n labels a set of density matrices in which 

Po = P{{o,o, - ,0} {o,o,---,o}} is the system reduced density operator (RDO), and all others 

are considered as auxiliary density operators (ADOs). Although the RDO is the most 
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important operator, the ADOs contain corrections to the system-bath interaction, arising 
from the non-equihbrium treatment of the bath. K is the truncation level for the correlation 
function (Matsuraba frequency and constant c^) and the cutoff for the tier of ADOs was 
set at A^. The scaled approach guarantees that all elements in the ADOs decay to zero for 
the upper levels in the hierarchy, while the Ishizaki-Tanimura truncating scheme decreases 
the truncation error. For a detailed derivation of this approach we refer the reader to |15] . 
We make use of the same parameters as [15j, and we set the truncation levels K = {] and 
cutoff tier of ADOs Mc = 4. The reorganization energy and Drude decay constant are 
A,- = A = 35 cm^-*^ and 7^"*^ = 7^-*^ = 50 fs. 



By numerically integrating the differential equation Eq. 12 using Mathematica, we cal- 
culated the density matrix of each time step during the evolution for 2500fs with a time 
step of 2 fs. We performed simulations with two different initial states: site 1 initially exited 
and site 6 initially excited. The time series of the system density matrix so obtained is the 
data from which we calculate the entanglement between various different parts of the FMO 
complex. Before describing the results of those calculations, we first describe the method 
by which they were obtained. 



III. ENTANGLEMENT ANALYSIS 



The FMO complex, considered as an assembly of seven chromophores, is a multipartite 
quantum system. As such, useful information about quantum correlations is obtained by 
computing the bipartite entanglement across any of the cuts that divide the seven chro- 
mophores into two subsystems. Similarly if we take the state of any subsystem of the FMO 
complex we can compute the entanglement across any cut of the reduced state of that 
subsystem. 



A. Entanglement measures 

In the present paper we choose to compute the measures of entanglement defined in |54] . 
which are closely related to previous measures defined by Meyer and Wallach [48j, Bren- 
nen |55], Scott |56] and Yu and Song |57]. A good starting point for consideration of these 
entanglement measures is the Meyer- Wallach measure defined in |48) . Given p = 
where {ip) G C^®" with K^IV")!^ = 1, the Meyer- Wallach measure can be written as fol- 
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lows mi [55]. 

1=1 

where is the reduced density matrix of the i*^ qubit. 

This measure is based on the fact that the purity or mixedness of the reduced density 
matrices pi can be viewed as an indicator of entanglement of the i^'^ qubit with the other 
n — 1: Tr(p^) < 1 such that equality is obtained if and only if p is a pure state. We are free 
to measure entanglement between any pair of subsystems in the same way. We use the set 
of monotones defined in ^54j : 

Vs = (1 - Tr(p|)) (14) 

where 5* is a set of k qubits, so that 15*1 = k, and ps is the reduced density matrix over those 
k qubits. 



B. Monogamy of entanglement 

A fascinating property distinguishing entanglement from classical correlations is 
monogamy. Just as the simplest example of entanglement occurs for two qubits, the simplest 
example of monogamy occurs for three qubits. If, among three qubits ABC, the qubits A 
and B are maximally entangled, then qubit C cannot be entangled at all with qubits A and 
B. It is instructive to consider this from the point of view of the entanglement measures (Eq. 



13). These measures are based on subsystem purity - if qubits ABC are in a pure state and 
A and B are maximally entangled then the reduced state of qubits AB is pure, hence so is 
the reduced state of qubit C, and hence qubit C is unentangled with qubits A and B. In fact, 
this property extends for three qubits to the case where the entanglement is not maximal. 
The monogamy constraint is expressed in terms of the tangles measuring the entanglement 
of qubit A with a subsystem B: 

TA\B = 2(1 - Trp^) = r,A. (15) 



In terms of the measures (Eq. 15) we obtain: 

taib + ta\c < ta\bc (16) 
This property of three qubit states was shown in [58], and the result for n qubits was proved 



m 



Tm\i < r{m\l, . . .m - l,m + 1 . . .n). (17) 
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These imply corresponding relations among the measures rjs that are equal to tangles of one 
qubit S with the others. 

In the context of models of exciton transport that are restricted to the single exciton 
subspace it is worth recalling that, in the case of pure states of three qubits, it is exactly states 
that are superpositions of Hamming weight one basis states that saturate the monogamy 
bound j5^. In fact, for both pure and mixed state cases it has been shown that generalized 
W states and mixtures of generalized W states with |0) (0| (which corresponds exactly to the 
case of interest for the single exciton subspace in the models we consider here) saturate the 
monogamy bounds [49j. We may therefore obtain the entanglement of each chromophore 
with the rest using the sum of the pairwise entanglements. The other entanglements, of 
pairs with pairs, and so on, may not be obtained from monogamy properties of qubit W- 
class states. 

It is natural to ask whether monogamy holds beyond restrictions on the entanglement 
of single qubits to relationships between the entanglement of higher dimensional systems. 
Unfortunately, this is not the case |60], as it can already be shown that states of qutrits 



violate the analogous relation to (Eq. 16). However, whether there are other inequalities 
among the full set of measures 7]s is not currently known. 



C. Convex Roof Extension of Entanglement Monotones 



Given a density matrix p and its set of ensemble representations 



(18) 



any entanglement monotone on pure states can be generalized to a monotone on 

mixed states, E{p), defined by 



(19) 



which is also an entanglement monotone. Given a density matrix p = YlPil'^di'^il^ define 

i 

\(|)^)Vq^ = J2^iJ\^^)VPJ, (20) 
j 

where the Ui/s are elements of a unitary matrix. It can then be shown that p = ^ q'j|</>i) {(f)i\- 

i 

Since density matrices are hermitian they are always diagonalizable. We can therefore 
write p = VAV''; this matrix product can equivalently be written as the summation p = 
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Xi\vi){vi\, where the Aj's are the eigenvalues of A and the |fi)'s are the basis-independent 

i 

orthonormal kets corresponding to the columns of V. This is called the spectral ensemble 
of p. It is also useful to define $ = VA^^'^, so that = p. This object $ contains 

all the information contained in a particular ensemble, and similar objects \E'\1'^ = p also 
correspond to ensembles. In fact, the unitary transform given in terms of a summation above 
corresponds to the matrix transformation $f/, where U is unitary. If we define \1/ = $f/ for 
some unitary matrix U, then ^^"^ = p. It can further be shown that the space of ensemble 
representations of p is isomorphic to the unitary group [61 j. Hence optimization over the 
space of ensembles can be reduced to an optimization problem over the unitary group. 

D. The Cayley Map 

The Cayley map is a self-inverse map from the algebra u{N) to the group U{N). The 
Cayley map is a map between a number of Lie algebras and their respective groups. It was 
introduced as a map from so{N) to SO{N) |62]. The Cayley map is defined by 

Cay(a) = A = {I - a) {I + a)^^ (21) 

where a is an element of the algebra being considered, and A is an element of the group. 
Likewise, we have 

CayA = a = {I-A){I + Ay^ (22) 

In the case of the unitary group, the Cayley map is a bijection between u{N) and the set 
U{N) — <f, where is the set of "exceptional elements." <f is the set of all elements A such 
that / + y4 is singular, and can be characterized as the set of all elements A with at least one 
eigenvalue —1. The exceptional elements on 5*0(3) are the reflections. For all such elements 
E, I + E has a eigenvalue, and is not invertible, so the Cayley map is not defined on these 
elements; however, this will not hinder our attempts to minimize 1] over U {N). Since we are 
performing numerical optimization, we only care that we can get arbitrarily close to a given 
local optimum. The closure of the image of the Cayley map on u{N) is all of U{N), so we 
will still be able to identify minima located at exceptional points. 

Because u{N) is easily parametrized by A^^ parameters, we can therefore parametrize 
U{N) by A^^ parameters via the Cayley map. Given a set of A^^ parameters {pi, . . . ,Pn'^}, 
the corresponding element of U{N) is then: 



A = Cay(a(pi,. . . ,Pn^)) 



(23) 
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where a is the element of u{N) given by the parameters pi under a standard parametrization. 
In the current work we use the basis of tensor products of Pauh matrices for the algebra 
su{N). The virtue of the Cayley map is that it gives us an easily understood and easily 
implemented way to parametrize U{N). The Cayley map thus provides somewhat simpler 
parameterization than that used in prior work on the convex roof optimization in |63]. Com- 
parison of the performance of our method with the simulated annealing approach described 
in Appendix B of [64J shows a substantial advantage to parameterization by the Cayley 
map combined with steepest descent. We leave detailed comparison of our method with 
that of [63j, and the evaluation of other optimization techniques beyond steepest descent, 
to future work. 

IV. RESULTS AND DISCUSSION 

There are 63 distinct bipartitions of the 7 chromophores of FMO. Ideally, one would com- 
pute all of these measures to obtain a complete picture of the correlations present among 
subsystems. In practice, one may use the saturation of the monogamy bounds to determine 
7 exactly, and the remaining 56 require computation of the convex roof. Instead one may 
take subsystems and compute the entanglement across bipartitions of the subsystems. For 
example, by computing the entanglement between all pairs of chromophores. However, as 
Table [T] shows, this leads to a large number of subsystems, and a large number of biparti- 
tions for each subsystem. Many of these measures are determined by the monogamy bounds 
in the single exciton subspace, but not all. In this section, we compute a number of en- 
tanglement measures for two, three, four and five qubit subsystems. Our approach follows 
both that of [42], in which pairwise entanglments were computed, and that of |46] in which 
the logarithmic negativity for several partitions of the full seven chromophore system were 
computed. The convex roof optimization described above is necessary, as the results of |42| 
cannot be extended beyond pairs without it as the optimal ensemble is not known analyti- 
cally beyond two qubits, and the results of j46j cannot be extended to measures other than 
logarithmic negativity because other measures do not share the computational tractability 
of the logarithmic negativity. 

A. Two chromophore subsystems 

The pairwise entanglements are a natural starting point because they may be computed 
exactly. For the case in which site 1 was initially excited the coherent oscillations of pop- 
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31 
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25 of 31 partitions not determined by monogamy 


7 


1 


63 


63 


7 


Seven determined from pairwise entanglements 



Table I: Subsystems and bipartite cuts relevant to the FMO system. One may take a subsystem 
reduced density matrix of any m < 7 and consider all the bipartite cuts of each subsystem. This 
leads to a combinatoric explosion of different bipartite measures. Evidently it would be simpler to 
consider all cuts of the total system, however, the cost of the convex roof optimization prohibits 
this at the present time. 



ulation occur mainly between sites 1 and 2 before the energy is transferred to sites 3 and 
4 [m [T5] . As a result of these coherent oscillations there is large pairwise entanglement 
between site 1|2 |l2]. In the work of |12], for times <900 (500) fs at 77K (300K) these 
measures are ordered: 

1|2 > 1|3 > 1|5 > 3|4 (24) 

For the system of pE] initialized with a single exciton at site 6 the entanglements 4|5, 4|7, 
5|6, 3|4 are computed. For times < 100 fs these are ordered 

5|6 > 4|5 > 4|7 > 3|4. (25) 

In Figure[T]we plot the entanglement evolution of the FMO complex when site 1 is initially 
excited at T = 77K. Figure [T] shows all 21 pairwise concurrences computed by the convex 
roof. For entanglements 1|2 and 1|3 we also plot the exact concurrence - the agreement is 
good enough that the difference between the convex roof and the exact calculation is not 
visible. In Figure |2] we plot the same data when site 6 is initially excited at T = 77K. 
For entanglements 5|6 and 4|5 we also plot the exact concurrence - again the agreement is 
good enough that the difference between the convex roof and the exact calculation is not 
visible. These show the ordering 1|2 > 1|3 > 1|5 as the significant entanglements for site 
one initially excited and 5|6 > 4|5 as the significant entanglements for site 6 initially excited 
consistent with the results of ||42jj. Because the monogamy bound is saturated in the single 



14 




Time / fs 

Figure 1: Entanglement evolution in the FMO complex when site 1 is initially excited at T = 77K. 
This Figure shows all 21 pairwise concurrences computed by the convex roof. For entanglements 
1|2 and 1|3 we also plot the exact concurrence - the agreement is good enough that the difference 
between the convex roof and thee exact calculation is not visible. Because the monogamy bound is 
saturated in the single exciton manifold, these 21 measures determine the entanglement of any one 
chromophore with any subset of the others. 

exciton manifold, these 21 measures determine the entanglement of any one chromophore 
with any subset of the others. 

These results on two chromophore subsystems help us identify a pathway involving sites 
1234 as significant for exciton transport when site 1 is initially excited, and a pathway 
involving sites 6543 as significant for exciton transport when site 6 is initially excited. This 
is consistent with prior results on pairwise entanglement |39t 112]. These results also validate 
our convex roof computations, at least for the case of two chromophore systems. It is perhaps 
unsuprising that the convex roof optimization performs well in that setting and so we now 
turn our attention to larger subsystems. 
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Figure 2: Entanglement evolution in the FMO complex when site 6 is initially excited at T = 77K. 
This Figure shows all 21 pairwise entanglements computed by the convex roof. For entanglements 
5|6 and 4|5 we also plot the exact concurrence - the agreement is good enough that the difference 
between the convex roof and thee exact calculation is not visible. Because the monogamy bound is 
saturated in the single exciton manifold, these 21 measures determine the entanglement of any one 
chromophore with any subset of the others. 



B. Three chromophore subsystems 

Figure |3] shows results beyond the pairwise entanglements. We compute the entanglement 
among the triples 134, 234, 123, 124 both from the monogamy bound and by the convex roof 
procedure. These results show both the utility of the saturation of the monogamy bound 
and the performance of the convex roof optimization. We see that the convex roof performs 
well for three qubits, closely matching the monogamy bound. 

Figure |4] shows the triple entanglement evolution among sites 123 in both the isolated 
system and the system with environmental coupling. The left side of the image shows the 
entanglement evolution for the isolated system, while the right side is the entanglement evo- 
lution under scaled HEOM approach. For the isolated system, the oscillations in population 
and entanglement will last forever. By comparing it with the open system case, it is obvious 
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Figure 3: Monogamy bound and convex roof computation of entanglements 1|34, 2|34, 12|3 and 
12|4. Particularly in the first 200 fs the convex roof closely matches the monogamy bound (which is 
equal to the corresponding tangle due to saturation of the monogamy bounds in the single exciton 
subspace). These calculations enable us to compute the entanglement among the significant triplets 
in FMO during transport. They also serve to validate the convex roof code which we shall use to 
compute entanglements not determined from the pairwise concurrences via monogamy. 



that the environment has the effect of eliminating the coherent oscillations characteristic of 
closed system quantum dynamics. Both isolated and system with environment case hit the 
maximum and minimum values at the same time during the evolution, which shows that 
the oscillations in the open system case are indeed the remnants of the coherent behavior 
in the closed system case. The entanglement evolution is not as smooth as ref. [42j, be- 
cause the simulation data has been sampled every lOfs in order to perform the entanglement 
calculations. 

Figjljo shows the entanglement of subsystem 123 across partition 1|23. The pairwise 
entanglement between site 1|2 and 1|3 and the monogamy bound is also shown. Because 
the monogamy bound is saturated for states in the zero and single exciton subspace 
these calculations show how well the convex roof optimization is performing. The time 
series of the triplet entanglement reflects the coherent oscillation of the population and 
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Figure 4: Entanglement evolution of FMO complex when site 1 is initially excited at cryogenic 
temperature T = 77K. The triplet site entanglement among site 1, 2 and 3 and also the pairwise 
site entanglement between any two of site 1, 2 and 3 are plotted. The left panel shows the dynamics 
of the entanglement for the system alone while the right considers the effect of the environment 

the time over which these oscillations last is the same as that in the population evolution 
which is around 650fs. The triplet entanglement evolution is comparable to the pairwise 
entanglement evolution between site 1|2. In the first a few beating region (t < 200fs), the 
triplet entanglement is almost the same as that of pairwise site 1|2. Beyond 200fs, the triplet 
entanglement becomes slightly larger than the pairwise entanglement site 1|2, indicating, via 
monogamy, that sites one and three have become entangled at this time. 

Fig. |4|l, shows the entanglement of the triplet 123 across the partition 2|13. This time 
series is similar to that of 1|23. Another interesting phenomena is the pairwise entanglement 
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between site 2|3, which also shows coherent oscillations. Although the amplitude of the 
pairwise entanglement is much smaller compared with the entanglement between site 1|2, 
they share the same frequency and hit the maximum and minimum value simultaneously. 

Fig. |4|[ shows the entanglement of the triplet 123 across the partition 3|12, which much 
smaller that the triplet entanglement 1|23 and 2|13 and does not show significant coher- 
ent oscillations. These entanglements are simply sums of pairwise entanglements, by the 
monogamy bound in the one-exciton subspace. For this case, in which site 1 is initially ex- 
cited, the dominant pairwise entanglement is 1|2, which is consistent with the other results 
in the literature jl2lli5| H6|. 

As a result, we conclude that in this pathway: during the coherent evolution period (first 
200fs), sites 3 and 4 are competing with each other to be entangled with site 12. However, 
when the coherent evolution disappears, the entanglement between site 3 and 4 becomes 
dominant. 

In order to check the thermal effect for the entanglement evolution, we plotted the entan- 
glement evolution under room temperature (T = 300T) for both site 1 and site 6 initially 
excited. By comparing with the evolution at T = 77K, the coherent oscillations were re- 
duced from 4 to 2 oscilations and the length of coherent oscillations was also reduced from 
650fs to 400fs. The maximum entanglement during the evolution was also reduced due to 
the increase in temperature. For example, the maximum triple entanglement of site 1|23 is 
0.85 at 77K while that is around 0.73 when T = 300K. In addition the entanglement goes 
to the equilibrium state much faster at 300K than at T = 77K. It takes around 7ps for the 
system to arrive at the equilibrium state at T = 77K, while at T = 300K this takes around 
1.5ps. The results at 300i^ are shown in Figure |5j 
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Figure 5: Time evolution of entanglement for multiple sites at T = 300K). The upper panel shows 
the entanglement when site 1 initially excited. Both triplet and pairwise sites entanglement among 
site 1, 2 and 3 are plotted. For the lower panel, site 6 is initially excited. We are focus on the 
quadruplet and pairwise entanglement among site 4 — 7. 

C. Four qubit subsystems 

Once again, for four qubit subsystems we may use the monogamy bounds to evaluate the 
performance of our convex roof calculations. 

In Figure |6] we evaluate the performance of our convex roof optimization using the 
monogamy bounds. As one can see, the agreement is less good than for two and three 
chromophore systems, but is significantly better in the case where the entanglements are 
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Figure 6: Entanglements in a four qubit system and monogamy bounds when site 1 is initially 
excited at temperature T = 77K. The entanglements 4|123, 7|123 and 7|456 are shown here - both 
computed exactly using the monogamy bound and determined by the convex roof. We see a larger 
variation in performance of the convex roof optimization here, with better agreement for 7|456 and 
7|123 than for 4|123. 

rather small. 

Next we examine the different roles of sites 3 and 4 in the pathway identified above for 
the case where site 1 is initially excited. It is known that the destination of this pathway is 
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Figure 7: Time evolution of entanglement related to site 3 and 4 for site 1 initially excited under 
T = 77K. Biparticle (site 3|4) triplet-particle (3|12; 4|12) and quadruplet-particle ( 4|123 ) are all 
considered here. 

the pair of sites 34. However, the detailed roles of these two sites during the entanglement 
evolution is still not clear. Figure [7] shows the entanglement evolution of the quadruple 1234 
across partition 4|123. It is clear that the pairwise entanglement 3|4 evolves in lockstep 
with 4 1 123 after 200fs. The entanglement of 3|12 and 4|12 are also evolving comparably 
after 200fs. Within the first 200 fs we see coherent oscillations in which 3|12 and 4|12 are 
in antiphase, but where 4|123 is in phase with 4|12. This behavior is suggestive of an initial 
period (the first 200 fs) in which the entanglement of chromophore 4 with 123 is fixed by its 
entanglement with chromophores 12, and then a long - time behavior in which chromophore 
4 is entangled with chromophore 3. This is consistent with a picture of energy transport in 
which a delocalized exciton passes from chromophores 12 to chromophores 34 - eventually 
landing at chromophore 3. 

In Figure |8] we compute the entanglement between pairs of chromophores 12 and 34 for 
the case where site 1 is initially excited. These entanglements are not bound by monogamy. 
Comparison of this figure with Figure [7] is instructive, as we see that the entanglement 
between the pairs of chromophores 12 and 34 is decreasing after the first 200 fs - following 
the falling entanglement of the pair 1|3. This makes sense in a picture of transport in which 
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Figure 8: Entanglements in a four qubit system when site 1 is initially excited at temperature 
T = 77K. The entanglement 12|34 computed via the convex roof procedure is shown here, together 
with the significant pairwise entanglements (computed exactly) among the quadruple. Although 
monogamy bounds do not apply in this case (as neither of the subsystems 12 or 34 is a qubit), we 
see that the entanglement 12|34 evolves similarly to the 1|3 entanglement. 

12 are the chromophores receiving the exciton when it is injected and 34 receive the exciton 
before it passes to the reaction center. 

For the case in which site 6 is initially excited, we compute entanglement among significant 
quadruplets of sites. Fig. |9]shows the entanglement evolution of the quadruplet 4567 for both 
the isolated and the open system case. Similar to the case where site 1 is initially excited, 
the entanglement displays coherent oscillations. Furthermore, these oscillations persist as 
long as the oscillations in the population. However, the timescale over which these coherent 
oscillations last is only 300fs, much shorter than the case where site 1 is initially excited. 
This phenomena is consistent with the population evolution [Til [15]. The most significant 
pairwise entanglement is 5|6, for which the maximum value can go as high as 0.8. The 
second most important pairs are sites 4|5 and 4|6, which have the maximum amplitude 
around 0.4. On the other hand, the coherent oscillations for all three pairs share the same 
frequency and evolution trend besides the 1st beating. The quadruplet entanglement 6|457 
and 5 1 467 have similar amplitude and evolution process, while the quadruplet 3|567 and 
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Figure 9: Entanglement evolution of FMO complex when site 6 is initially excited at temperature 
T = 77K. The quadruplet sites entanglement for site 4—7 and the bi-site entanglement between 
each two sites among them are plotted. The left panel shows the isolated situation and the right 
panel considers the dynamics with environment. 
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4|567 are much smaller compared with the above two. During the first two beatings, the 
triplet entanglement is almost the same as the pairwise entanglement. This shows again 
that the triplet entanglement under this initial condition is still dominated by the pairwise 
entanglement - consistent with the saturation of monogamy bounds in this case. 



D. Five qubit subsystems 
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Figure 10: Entanglement evolution of the quintuplet sites ( 5|1234, 6|1234 and 7|1234 ) and the 
corresponding monogamy bounds in the FMO complex under the cryogenic temperature T = 77K. 
Site 1 is initially excited. Site 1, 2, 3 and 4 are sites evolved in the population pathway under this 
initial condition. 



Fig. [To] shows the quintuplet entanglement evolution among site 5, 6 and 7 with quadru- 
plet sites 1234. All three combinations 5|1234, 6|1234 and 7|1234 have very small entan- 
glement during the time evolution. For the case in which site 1 is initially excited, the 
entanglement evolution only happens among sites in the pathway, which is site 1, 2, 3 and 4. 
We also plotted the monogamy bounds in Fig. [10} which indicates the monogamy bounds 
are always smaller than the corresponding quintuplet entanglements - this shows that, un- 
suprizingly, the convex roof optimization is not performing as well in the five qubit case as 
it does for three and four qubits. 

For the case in which site 1 is initially excited, entanglement is only significant within the 
sites in the pathway. We would like to know if this is also the case when site 6 is initially 



excited. Figure 11 shows the quintuplet entanglement for sites 1 and 2 with sites 4567. 
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Figure 11: Time evolution of entanglement among five-site 1|4567 and 2 1 4567 in the FMO complex 
for site 6 initially excited situation. The system is under cryogenic temperature 77K 

The maximum entanglement for those 2 is around 0.25, which is much smaller compared 
with that among quadruplet sites 4567. This is consistent with the idea that entanglement 
is concentrated among the sites evolved in a specific pathway, with different pathways for 
different initial conditions. 

As for the site 1 initially excited case, we also examined the roles of sites 3 and 4 in the 



case when site 6 is initially excited (Fig. 12). Just as in the case where sit 1 was initially 
excited (Figure [T]) we see a 200fs period with coherent oscillations in the entanglement in 
which the entanglement of 3 with the rest and 4 with the rest are in antiphase. This is 
followed by a later period in which sites 3 and 4 become entangled and the entanglement 
of 3 with 4567 is dominated by the entanglement of 3 and 4. As a result, the dominant 
pairwise entanglement changes from site 5|6 to pair 3|4 during the transport of the exciton 
from the injection site at site 6 to the final state in which it is concentrated on sites 3 and 
4. 
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Figure 12: Similar to Figure [Tj the entanglement evolution related to site 3 and 4 are investigated 
when site 6 is initially excited. The entanglement of site 3 and 4 with the other three sites 5, 6 and 
7 is considered 

E. Beyond the single exciton manifold 

In addition to measuring entanglement in the one-exciton subspace, we conducted a 
number of tests where we manually reinserted the ground state density matrix po = 
|0000000)(0000000| and the two-exciton density matrix pa = |0000011)(0000011| in order 
to try to determine how the entanglement would be affected. In the first test, we inserted 
the ground state po on its own, yielding the following expression for the density matrix 
(where pi is the density matrix for the single-exciton subspace): 



In our second test, we added in the two-exciton subspace alone, without the ground state: 

P = '-^P^^ (27) 



We then inserted po and p2 as follows: 

P = 1 , 1^,12 , 1^,14 /o (28) 
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Figure 13: A comparison of the effects of adding in the two-exciton subspace for different values of 
|a^|. The entanglement between sites one and two is plotted for the density matrix in equation 
with la^l = .5 .1 .01. 
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When we added in both the vacuum state and the two-exciton subspace |0000011)(0000011| 
and varied a, we found that for values as small as \a\'^ = .01, the entanglement completely 
disappeared. We then experimented with adding in both the ground state and an exponen- 
tially decaying two-exciton subspace, p2 = e~'^*|0000011)(0000011|, and, as expected, as e~'^* 
goes to zero, we recover some entanglement between sites 1 and 2, although the magnitude 



is still diminished by the presence of the vacuum state 29 In order to get a sense of how 
quickly the entanglement recovers, we calculated the concurrence for the density matrix 
in equation [28| which includes the ground state 1 0000000) (0000000 1 and the two-exciton 
subspace |0000011)(0000011| scaled by a factor 7 e [0, 1]: 

Po + + 7l«l^/2p2 



1 + + 7|a|V2 



(29) 



The results are plotted in Figs. 14 and 15 
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Figure 14: A comparison of the effects of adding in the two-exciton subspace 

7|ap|0000011) (OOOOOllj for different values of 7, with |ap = 0.5. The entanglement between 



sites one and two is plotted for the density matrix in 28 with p2 = |0000011)(0000011 



V. CONCLUSIONS 



In summary, we used the direct computation of the convex roof optimization approach to 
simulate the entanglement evolution in the FMO complex via the scaled HEOM approach. 
Our simulation results indicate the following results: 

Pairwise entanglement plays the dominant role in entanglement of the FMO complex. 
Because of the saturation of the monogamy bounds the entanglement of any chromophore 
with any subset of the other chromophores is completely determined by the set of pairwise 
entanglements. For the simulations in which site 1 is initially excited, the dominant pair is 
site 1 and 2, while in the cases where 6 is initially excited site 5 and 6 are most entangled. 
This indicates that entanglement is dominant in the early stages of exciton transport, when 
the exciton is initially delocalized away from the injection site. In addition we observe that 
the entanglement mainly happens among the sites involved in the pathway. For the site 
1 initially excited case, the entanglement of site 5, 6 and 7 is almost zero. For the site 6 
initially excited situation, there is seldom entanglement for site 1 and 2. 
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Figure 15: A curve showing how the amphtude of the entanglement between sites 1 and 2 at 21 



fs varies as a function of 7, with the density matrix from 28 with p2 = |0000011)(0000011| and 
lap = .5. 



Although the final state is the same for both initial conditions, the role of site 3 and site 
4 during the time evolution is different. For the initial condition where site 1 is excited, the 
entanglement is transferred to site 3 and then from site 3 to site 4. While for the site 6 
initially excited case, sites 4 and 5 first become entangled with site 6 and then sites 3 and 
4 become entangled. This is due to the fact that site 3 has strong coupling with site 1 and 
2, while site 4 is coupled more strongly to sites 5, 6 and 7. 

The initial condition plays an important role in the entanglement evolution, the entan- 
glement decays faster for the cases where site 6 is initially excited compared with cases 
where the site 1 is initially excited. Increasing the temperature unsuprisingly reduces the 
amplitude of the entanglement and also decreases the time for the system goes to thermal 
equilibrium. 
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VI. FUTURE DIRECTIONS 

Most entanglement measures computed previously for FMO were chosen on the basis of 
ease of calculation. The negativity and logarithmic negativity are straightforward to com- 
pute for all states [361 146] . The global and bipartite relative entropy of entanglement can be 
made straightforward to compute by restriction to the single exciton subspace |^ 151). The 
bipartite concurrence and tangles can be computed easily for pairs of chromophores |42t H7]. 
In all cases the chosen measures of simplifications thereof enable one to avoid computing the 
convex roof over different ensembles representing a mixed state. In this paper we explored 
the difficulty of such calculations, and find that measures that yield the bipartite entangle- 
ment across cuts of 3,4, and 5 qubit subsystems may be computed with modest effort. We 
computed monogamy bounds, using the saturation of these bounds in the single exciton sub- 
space to evaluate the performance of our convex roof procedure. This technique enables 
us to extend the set of measures that have been computed for FMO, and also shows that the 
computation of entanglement for this system is not restricted by the difficulty of the convex 
roof procedure. This procedure could also be used, with no increase in computationla cost, 
to analyze entanglement in multiexcitonic models. 
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